rm(list = ls())
###

###
library(tidyverse) # 1.3.1
library(stargazer) # 5.2.2
library(scales) # 1.1.1
library(modelsummary) # 0.7.0
library(vtable) # 1.3.1
library(gtsummary) # 1.4.1
library(gridExtra) # 2.3
library(jtools) # 2.1.3

## Read in ACS Data
acs19 <- haven::read_dta("acs_2019.dta")

## Code variables
acs19 %>% 
  filter(age >= 18,
         ftotinc >= 0) %>%
  select(sex, age, region, ftotinc, educd, race, rachsing, perwt) %>%
  mutate(gender = case_when(sex == 1 ~ "Male",
                            sex == 2 ~ "Women"),
         age = case_when(age >= 18 & age <= 29 ~ "18 - 29",
                         age >= 30 & age <= 41 ~ "30 - 41",
                         age >= 42 & age <= 54 ~ "42 - 54",
                         age >= 55 ~ "55 or more",
                         TRUE ~ "other"),
         region = case_when(region >= 11 & region <= 13 ~ "Northeast",
                            region >= 21 & region <= 23 ~ "Midwest",
                            region >= 31 & region <= 34 ~ "South",
                            region >= 41 & region <= 43 ~ "West"),
         educ = case_when(educd == 00 | educd == 02 |
                            educd >= 003 & educd <= 061 ~ "Less than HS",
                          educd >= 062 & educd <= 064 ~ "HS",
                          educd >= 065 & educd <= 080 |
                            educd == 090 | educd == 100 | 
                            educd >= 110 & educd <= 113 ~ "Some college",
                          educd >= 081 & educd <= 083 ~ "AA",
                          educd == 101 ~ "BA",
                          educd >= 114 & educd <= 116 ~ "Postgraduate",
                          educd == 999 | educd == 001 ~ "NA",
                          TRUE ~ "NA"),
         education = case_when(educd == 000 | educd == 002 |
                                 educd >= 003 & educd <= 061 ~ 1,
                               educd >= 062 & educd <= 064 ~ 2,
                               educd >= 065 & educd <= 080 |
                                 educd == 090 | educd == 100 | 
                                 educd >= 110 & educd <= 113 ~ 3,
                               educd >= 081 & educd <= 083 ~ 4,
                               educd == 101 ~ 5,
                               educd >= 114 & educd <= 116 ~ 6,
                               educd == 999 | educd == 001 ~ 999),
         income = case_when(ftotinc >= 0 & ftotinc <= 20000 ~ "20k or less",
                            ftotinc > 20000 & ftotinc <= 59999 ~ "20 - 59k",
                            ftotinc >= 60000 & ftotinc <= 99999 ~ "60 - 99k",
                            ftotinc >= 100000 & ftotinc <= 9999997 ~ "100k or more",
                            ftotinc == 9999998 | ftotinc == 9999999 ~ "N/A"),
         rachsing2 = case_when(rachsing == 1 ~ "White",
                               rachsing == 2 ~ "Black",
                               rachsing == 3 ~ "Native American",
                               rachsing == 4 ~ "Asian",
                               rachsing == 5 ~ "Hispanic")) %>%
  rename(ethnicity  = rachsing2) %>%
  select(gender, educ,  age, region, income, ethnicity) -> descrip_acs19


# Add indicator
descrip_acs19 %>%
  mutate(study = rep(3, length = nrow(descrip_acs19)),
         gender = as.factor(gender),
         ethnicity = as.factor(ethnicity),
         educ = as.factor(educ),
         age = as.factor(age),
         region = as.factor(region),
         income = as.factor(income)) -> descrip_acs19


#############################################     
########## Load Survey ######################
#############################################     

######
data <- read.csv("perspective_taking_merged.csv")

###### Descriptives of merged survey
data <- data %>%
  filter(informed_consent == 1 & # Remove refused informed consent
           attention_check1 == 3 &
           age >= 18) %>% # Remove failed attention check #1
  mutate(ethnicity = factor(ethnicity),
         income = factor(case_when(
           income == 1 ~ "20k or less",
           income >= 2 & income <= 6 ~ "20 - 59k",
           income >= 7 & income <= 9 ~ "60 - 99k",
           income == 12 ~ "100k or more",
           TRUE ~ "N/A")),
         region = factor(case_when(
           region == 1 ~ "Northeast",
           region == 2 ~ "Midwest",
           region == 3 ~ "South",
           region == 4 ~ "West")),
         gender = factor(case_when(gender == 1 ~ "Male",
                            gender == 2 ~ "Women")),
         educ = factor(case_when(
           educ == 1 ~ "Less than HS", 
           educ == 2 ~ "HS",
           educ == 3 ~ "Some college",
           educ == 4 ~ "AA",
           educ == 5 ~ "BA",
           educ == 6 ~ "Postgraduate",
           TRUE ~ "N/A")),
         age = factor(case_when(
           age >= 18 & age <= 29 ~ "18 - 29",
           age >= 30 & age <= 41 ~ "30 - 41",
           age >= 42 & age <= 54 ~ "42 - 54",
           age >= 55 ~ "55 or more",
           TRUE ~ "N/A")),
         ethnicity = case_when(
           ethnicity == 1 ~ "White",
           ethnicity == 2 ~ "Hispanic",
           ethnicity == 3 ~ "Black",
           ethnicity == 4 ~ "Asian",
           ethnicity == 5 ~ "Native American",
           TRUE ~ "other"),
         pid_3level = factor(case_when(pid_7level == 1 ~ 1,
                                pid_7level == 2 ~ 1,
                                pid_7level == 3 ~ 1,
                                pid_7level == 4 ~ 2,
                                pid_7level == 5 ~ 3,
                                pid_7level == 6 ~ 3,
                                pid_7level == 7 ~ 3)),
         polid_3level = factor(case_when(pol_ideology == 1 ~ 1,
                                  pol_ideology == 2 ~ 1,
                                  pol_ideology == 3 ~ 1,
                                  pol_ideology == 4 ~ 2,
                                  pol_ideology == 5 ~ 3,
                                  pol_ideology == 6 ~ 3,
                                  pol_ideology == 7 ~ 3)),
         covid_ever = factor(ifelse(covid_contracted == 1 |
                               covid_self == 1, "Yes", "No")) ) 

## select only variables you want
data %>%
  select(study, gender, educ, age, region, ethnicity, income) %>%
  mutate(gender = as.factor(gender),
         educ = as.factor(educ),
         region = as.factor(region),
         ethnicity = as.factor(ethnicity),
         income = as.factor(income)) -> study_merge


#### Join
descrip_all <- rbind(study_merge, descrip_acs19) %>% 
  mutate(educ = factor(educ,
                       levels = c("Less than HS", "HS",
                                  "Some college","AA",
                                  "BA","Postgraduate")),
         income = factor(income,
                        levels = c("20k or less", "20 - 59k",
                                   "60 - 99k", "100k or more", 
                                   "N/A")),
         study = factor(case_when(study == 1 ~ "Study 1",
                                  study == 2 ~ "Study 2",
                                  study == 3 ~ "ACS 2019"),
                        levels = c("Study 1",
                                   "Study 2",
                                   "ACS 2019")) )


########################################
########## Table A1 ####################
########################################

sumtable(descrip_all, 
         group = "study",
         out = "viewer")


########################################
########## Table A2 ####################
########################################

#### Survey Exp 1
data %>%
  filter(study == 1) %>%
  mutate(treatment = factor(treatment,
                            levels = c('Control Low Inc',
                                       'Control Undoc',
                                       'Treat Low Inc',
                                       'Treat Undoc')),
         educ = factor(educ,
                       levels = c("Less than HS", "HS",
                                  "Some college","AA",
                                  "BA","Postgraduate")),
         income = factor(income,
                         levels = c("20k or less", "20 - 59k",
                                    "60 - 99k", "100k or more",
                                    "N/A")),
         pid_3level = factor(case_when(pid_3level == 1 ~ "Democrat",
                                       pid_3level == 2 ~ "Independent",
                                       pid_3level == 3 ~ "Republican"),
                             levels = c("Democrat", "Independent", "Republican")),
         polid_3level = factor(case_when(polid_3level == 1 ~ "Liberal",
                                         polid_3level == 2 ~ "Moderate",
                                         polid_3level == 3 ~ "Conservative"),
                               levels = c("Liberal", "Moderate", "Conservative")) ) %>%
  select(treatment, gender, educ, age,
         region, ethnicity, income,
         pid_3level, polid_3level) %>%
  sumtable(group = "treatment",
           out = "viewer")

########################################
########## Table A3 ####################
########################################

#### Survey Exp 2
data %>%
  filter(study == 2) %>%
  mutate(educ = factor(educ,
                       levels = c("Less than HS", "HS",
                                  "Some college","AA",
                                  "BA","Postgraduate")),
         income = factor(income,
                         levels = c("20k or less", "20 - 59k",
                                    "60 - 99k", "100k or more", 
                                    "N/A")),
         treatment = factor(treatment,
                     levels = c("Control Undoc",
                                "Treat Undoc",
                                "Cancer",
                                "Immigrant")),
  pid_3level = factor(case_when(pid_3level == 1 ~ "Democrat",
                                pid_3level == 2 ~ "Independent",
                                pid_3level == 3 ~ "Republican"),
                      levels = c("Democrat", "Independent", "Republican")),
  polid_3level = factor(case_when(polid_3level == 1 ~ "Liberal",
                                  polid_3level == 2 ~ "Moderate",
                                  polid_3level == 3 ~ "Conservative"),
                        levels = c("Liberal", "Moderate", "Conservative")) ) %>%
  select(treatment, gender, educ, age,
         region, ethnicity, income,
         pid_3level, polid_3level) %>%
  sumtable(group = "treatment",
           out = "viewer")


########################################
########## Table A4 ####################
########################################

# Study 1
data %>%
  filter(study == 1) %>%
  mutate(ethnicity = factor(ethnicity),
         pid_3level = factor(pid_3level),
         polid_3level = factor(polid_3level),
         educ = factor(educ,
                        levels = c("Less than HS", "HS",
                                   "Some college","AA",
                                   "BA","Postgraduate")),
         income = factor(income,
                         levels = c("20k or less", "20 - 59k",
                                    "60 - 99k", "100k or more", 
                                    "N/A")),
         pid_3level = factor(case_when(pid_3level == 1 ~ "Democrat",
                                       pid_3level == 2 ~ "Independent",
                                       pid_3level == 3 ~ "Republican"),
                             levels = c("Democrat", "Independent", "Republican")),
         polid_3level = factor(case_when(polid_3level == 1 ~ "Liberal",
                                         polid_3level == 2 ~ "Moderate",
                                         polid_3level == 3 ~ "Conservative"),
                               levels = c("Liberal", "Moderate", "Conservative")) ) %>%
  select(covid_ever, gender, age,
         educ, ethnicity, income, 
         pid_3level, polid_3level) %>%
  sumtable(group = "covid_ever",
           out = "viewer")

########################################
########## Table A5 ####################
########################################

# Study 2
data %>%
  filter(study == 2) %>%
  mutate(ethnicity = factor(ethnicity),
         pid_3level = factor(pid_3level),
         polid_3level = factor(polid_3level),
         educ = factor(educ,
                       levels = c("Less than HS", "HS",
                                  "Some college","AA",
                                  "BA","Postgraduate")),
         income = factor(income,
                         levels = c("20k or less", "20 - 59k",
                                    "60 - 99k", "100k or more", 
                                    "N/A")),
         pid_3level = factor(case_when(pid_3level == 1 ~ "Democrat",
                                       pid_3level == 2 ~ "Independent",
                                       pid_3level == 3 ~ "Republican"),
                             levels = c("Democrat", "Independent", "Republican")),
         polid_3level = factor(case_when(polid_3level == 1 ~ "Liberal",
                                         polid_3level == 2 ~ "Moderate",
                                         polid_3level == 3 ~ "Conservative"),
                               levels = c("Liberal", "Moderate", "Conservative")) ) %>%
  select(covid_ever, gender, age,
         educ, ethnicity, income, 
         pid_3level, polid_3level) %>%
  sumtable(group = "covid_ever", 
           out = "viewer")


